/* Copyright (c) Colorado School of Mines, 2010.*/
/* All rights reserved.                       */

/* SUMIGTOPO2D: $Revision: 1.5 $ ; $Date: 2006/11/07 22:58:42 $	*/

#include "su.h"
#include "segy.h"

/*********************** self documentation **********************/
char *sdoc[] = {
" 									",
" SUMIGTOPO2D - Kirchhoff Depth Migration of 2D postack/prestack data	",
"	     from the (variable topography) recording surface		",
" 									",
"    sumigtopo2d  infile=  outfile=  [parameters] 			",
"									",
" Required parameters:							",
" infile=stdin		file for input seismic traces			",
" outfile=stdout	file for common offset migration output  	",
" ttfile		file for input traveltime tables		",
"   The following 9 parameters describe traveltime tables:		",
" fzt 			first depth sample in traveltime table		",
" nzt 			number of depth samples in traveltime table	",
" dzt			depth interval in traveltime table		",
" fxt			first lateral sample in traveltime table	",
" nxt			number of lateral samples in traveltime table	",
" dxt			lateral interval in traveltime table		",
" fs 			x-coordinate of first source			",
" ns 			number of sources				",
" ds 			x-coordinate increment of sources		",
"									",
" fxi                   x-coordinate of the first input trace           ",
" dxi                   horizontal spacing of input data                ",
" nxi                   number of input trace locations in surface      ",
"									",
" Optional Parameters:							",
" dt= or from header (dt) 	time sampling interval of input data	",
" ft= or from header (ft) 	first time sample of input data		",
" dxm= or from header (d2) 	sampling interval of midpoints 		",
" surf=\"0,0;99999,0\"  Recording surface \"x1,z1;x2,z2;x3,z3;...\"     ",
" fzo=fzt               z-coordinate of first point in output trace 	",
" dzo=0.2*dzt		vertical spacing of output trace 		",
" nzo=5*(nzt-1)+1 	number of points in output trace		",	
" fxo=fxt               x-coordinate of first output trace 		",
" dxo=0.5*dxt		horizontal spacing of output trace 		",
" nxo=2*(nxt-1)+1  	number of output traces 			",	
" off0=0               	first offest in output 				",
" doff=99999		offset increment in output 			",
" noff=1       		number of offsets in output 			",	
" fmax=0.25/dt		frequency-highcut for input traces		",
" offmax=99999		maximum absolute offset allowed in migration 	",
" aperx=nxt*dxt/2  	migration lateral aperature 			",
" angmax=60		migration angle aperature from vertical 	",
" v0=1500(m/s)		reference velocity value at surface		",	
" dvz=0.0  		reference velocity vertical gradient		",
" ls=1	                flag for line source				",
" jpfile=stderr		job print file name 				",
" mtr=100  		print verbal information at every mtr traces	",
" ntr=100000		maximum number of input traces to be migrated	",
"									",
" Notes:								",
" 1. Traveltime tables were generated by program rayt2dtopo (or any 	",
"    other one that considers topography )on relatively coarse grids,	",
"    with dimension ns*nxt*nzt. In the migration process, traveltimes	",
"    are interpolated into shot/gephone positions and output grids.	",
" 2. Input seismic traces must be SU format and can be any type of 	",
"    gathers (common shot, common offset, common CDP, and so on).	", 
" 3. Migrated traces are output in CDP gathers if velocity analysis	",
"    is required, with dimension nxo*noff*nzo.  			", 
" 4. If the offset value of an input trace is not in the offset array 	",
"    of output, the nearest one in the array is chosen. 		",
" 5. Amplitudes are computed using the reference velocity profile, v(z),",
"    specified by the parameters v0= and dvz=.				",
" 6. Input traces must specify source and receiver positions via the header",
"    fields tr.sx and tr.gx. Offset is computed automatically.		",
"									",
NULL};
/*
 * Author:  Zhenyue Liu, 03/01/95,  Colorado School of Mines
 *
 *	    Trino Salinas, 07/01/96, Colorado School of Mines,
 *          extended the code to migrate data from a nonflat 
 *          recording surface.
 *
 * References :
 *
 * Bleistein, N., Cohen, J., and Hagin, F., 1987, Two and one-half
 *   dimensional Born inversion with arbitrary reference: Geophysics,
 *   52, 26-36.
 *
 * Liu,Z., 1993, A Kirchhoff approach to seismic modeling and 
 *   pre-stack depth migration: CWP Annual Report, CWP, Colorado
 *   School of Mines.
 *
 * Wiggins, J. W., 1984, Kirchhoff integral extrapolation and
 *   migration of nonplanar data: Geophysics, 49, 1239-
 *   1248.
 *
**************** end self doc ***********************************/
/* TYPEDEFS */
/* Geometry of the recording surface */
typedef struct SurfaceSegmentStruct {
        float x;        /* x coordinate of segment midpoint */
        float z;        /* z coordinate of segment midpoint */
        float s;        /* x component of unit-normal-vector */
        float c;        /* z component of unit-normal-vector */
} SurfaceSegment;
typedef struct SurfaceStruct {
        int ns;               /* number of surface segments */
        float ds;             /* segment length */
        SurfaceSegment *ss;   /* array[ns] of surface segments */
} Surface;

/* Prototypes of functions used internally */
  void resit(int nx,float fx,float dx,int nz,int nr,float dr,
	float **tb,float **t,float x0);
  void interpx(int nxt,float fxt,float dxt,int nx,float fx,float dx,
	int nzt,float **tt,float **t);
  void sum2(int nx,int nz,float a1,float a2,float **t1,float **t2,float **t);
  void timeb(int nr,int nz,float dr,float dz,float fz,float sz,float a,
        float v0,float **t,float **p,float **sig,float **ang);
  void mig2d(float *trace,int nt,float ft,float dt,
	float sx,float gx,float **mig,float aperx,
  	int nx,float fx,float dx,float nz,float fz,float dz,
	int ls,int mtmax,float dxm,float fmax,float angmax,
	float **tbs,float **tbg,float **pbs,
	float **pbg,float **angbs,float **angbg,int nr,float **tsum, 
	int nzt,float fzt,float dzt,int nxt,float fxt,float dxt,float *szif,
	float nangls,float nanglg,float **sigbs,float **sigbg);
  void decodeSurfaces(int *nrPtr,int **nxzPtr,float ***xPtr,float ***zPtr);
  int decodeSurface(char *string,int *nxzPtr,float **xPtr,float **zPtr);
  void makesurf(float dsmax,int nr,int *nu,float **xu,float **zu,
        Surface **r);
  void zcoorSurfaces(float fx,float dx,int nx,float fxs,float dxs,int nxs,
        Surface *srf,float *szif,float *sz,float *nangl);

/* segy trace */
segy tr, tro;

int
main (int argc, char **argv)
{
	int   nt,nzt,nxt,nzo,nxo,ns,noff,nr,is,io,ixo,izo;
	int   ls,ntr,jtr,ktr,mtr,mtmax,nsrf,*nxzsrf,nxi;
	off_t nseek;
	float ft,fzt,fxt,fzo,fxo,fs,off0,dt,dzt,dxt,dzo,dxo,ds,doff,dxm,
	      ext,ezt,ezo,exo,es,s,scal,**xsrf,**zsrf,fxi,dxi;
	float v0,v00,dvz,fmax,angmax,offmax,rmax,aperx,sx,gx,nangls,
	      nanglg;
	float ***mig,***ttab,***tb,***pb,***angb,**tsum,**tt,*szif,
              *nangl,*sz,**tbs,**tbg,**pbs,**pbg,**angbs,**angbg,
              ***sigb,**sigbs,**sigbg;
        Surface *srf;
	
	char *infile="stdin",*outfile="stdout",*ttfile,*jpfile;
	FILE *infp,*outfp,*ttfp,*jpfp;


	/* hook up getpar to handle the parameters */
	initargs(argc, argv);
	requestdoc(1);

	/* open input and output files	*/
	if( !getparstring("infile",&infile)) {
		infp = stdin;
	} else  
		if ((infp=fopen(infile,"r"))==NULL)
			err("cannot open infile=%s\n",infile);
	if( !getparstring("outfile",&outfile)) {
		outfp = stdout;
	} else {
		outfp = fopen(outfile,"w");
	}

	efseeko(infp,(off_t) 0, SEEK_CUR);
	efseeko(outfp,(off_t) 0,SEEK_CUR);
	if( !getparstring("ttfile",&ttfile))
		err("must specify ttfile!\n");
	if ((ttfp=fopen(ttfile,"r"))==NULL)
		err("cannot open ttfile=%s\n",ttfile);
	if( !getparstring("jpfile",&jpfile)) {
		jpfp = stderr;
	} else  
		jpfp = fopen(jpfile,"w");

	/* get information from the first header */
	if (!fgettr(infp,&tr)) err("can't get first trace");
	nt = tr.ns;
	if (!getparfloat("dt",&dt)) dt = tr.dt/1000000.0; 
	if (dt<0.0000001) err("dt must be positive!\n");
	if (!getparfloat("ft",&ft)) ft = (float)tr.delrt/1000.; 
 	if (!getparfloat("dxm",&dxm)) dxm = tr.d2;
	if  (dxm<0.0000001) err("dxm must be positive!\n");
	
	/* get traveltime tabel parameters	*/
	if (!getparint("nxt",&nxt)) err("must specify nxt!\n");
	if (!getparfloat("fxt",&fxt)) err("must specify fxt!\n");
	if (!getparfloat("dxt",&dxt)) err("must specify dxt!\n");
	if (!getparint("nzt",&nzt)) err("must specify nzt!\n");
	if (!getparfloat("fzt",&fzt)) err("must specify fzt!\n");
	if (!getparfloat("dzt",&dzt)) err("must specify dzt!\n");
	if (!getparint("ns",&ns)) err("must specify ns!\n");
	if (!getparfloat("fs",&fs)) err("must specify fs!\n");
	if (!getparfloat("ds",&ds)) err("must specify ds!\n");
	ext = fxt+(nxt-1)*dxt;
	ezt = fzt+(nzt-1)*dzt;
	es = fs+(ns-1)*ds;

        if (!getparint("nxi",&nxi)) err("must specify nxi!\n");
        if (!getparfloat("fxi",&fxi)) err("must specify fxi!\n");
        if (!getparfloat("dxi",&dxi)) err("must specify dxi!\n");

	/* optional parameters	*/
	if (!getparint("nxo",&nxo)) nxo = (nxt-1)*2+1;
	if (!getparfloat("fxo",&fxo)) fxo = fxt;
	if (!getparfloat("dxo",&dxo)) dxo = dxt*0.5;
	if (!getparint("nzo",&nzo)) nzo = (nzt-1)*5+1;
	if (!getparfloat("fzo",&fzo)) fzo = fzt;
	if (!getparfloat("dzo",&dzo)) dzo = dzt*0.2;
	exo = fxo+(nxo-1)*dxo;
	ezo = fzo+(nzo-1)*dzo;
	if(fxt>fxo || ext<exo || fzt>fzo || ezt<ezo) 
		err(" migration output range is out of traveltime table!\n");

	if (!getparfloat("v0",&v0)) v0 = 1500;
	if (!getparfloat("dvz",&dvz)) dvz = 0;
	if (!getparfloat("angmax",&angmax)) angmax = 60.;
	if  (angmax<0.00001) err("angmax must be positive!\n");
	mtmax = 2*dxm*sin(angmax*PI/180.)/(v0*dt);
	if(mtmax<1) mtmax = 1;
	if (!getparfloat("aperx",&aperx)) aperx = 0.5*nxt*dxt;
	if (!getparfloat("offmax",&offmax)) offmax = 3000;
	if (!getparfloat("fmax",&fmax)) fmax = 0.25/dt;
	if (!getparint("noff",&noff))	noff = 1;
	if (!getparfloat("off0",&off0)) off0 = 0.;
	if (!getparfloat("doff",&doff)) doff = 99999;

	if (!getparint("ls",&ls))	ls = 1;
	if (!getparint("ntr",&ntr))	ntr = 100000;
	if (!getparint("mtr",&mtr))	mtr = 100;

        /* Estimation of the recording surface  */
        decodeSurfaces(&nsrf,&nxzsrf,&xsrf,&zsrf);
        makesurf(dxm*0.025,nsrf,nxzsrf,xsrf,zsrf,&srf);

	fprintf(jpfp,"\n");
	fprintf(jpfp," Migration parameters\n");
	fprintf(jpfp," ================\n");
	fprintf(jpfp," infile=%s \n",infile);
	fprintf(jpfp," outfile=%s \n",outfile);
	fprintf(jpfp," ttfile=%s \n",ttfile);
	fprintf(jpfp," \n");
	fprintf(jpfp," nzt=%d fzt=%g dzt=%g\n",nzt,fzt,dzt);
	fprintf(jpfp," nxt=%d fxt=%g dxt=%g\n",nxt,fxt,dxt);
 	fprintf(jpfp," ns=%d fs=%g ds=%g\n",ns,fs,ds);
	fprintf(jpfp," \n");
        fprintf(jpfp," nxi=%d fxi=%g dxi=%g\n",nxi,fxi,dxi);
        fprintf(jpfp," \n");
	fprintf(jpfp," nzo=%d fzo=%g dzo=%g\n",nzo,fzo,dzo);
	fprintf(jpfp," nxo=%d fxo=%g dxo=%g\n",nxo,fxo,dxo);
	fprintf(jpfp," \n");

	/* compute the z coordinate of the recording surface  */
        szif = ealloc1float(nxi);
        sz = ealloc1float(ns);
        nangl = ealloc1float(nxi);
        zcoorSurfaces(fxi,dxi,nxi,fs,ds,ns,srf,szif,sz,nangl);
	
	/* compute reference traveltime and slowness  */
	rmax = MAX(es-fxt,ext-fs);
	rmax = MIN(rmax,0.5*offmax+aperx);
	nr = 2+(int)(rmax/dxo);
        tb = ealloc3float(nzt,nr,ns);
        pb = ealloc3float(nzt,nr,ns);
        sigb = ealloc3float(nzt,nr,ns);
        angb = ealloc3float(nzt,nr,ns);
        for (is=0;is<ns;is++){
            v00 = v0 + sz[is]*dvz;
            timeb(nr,nzt,dxo,dzt,fzt,sz[is],dvz,v00,tb[is],pb[is],sigb[is],
                angb[is]);
        }

	fprintf(jpfp," nt=%d ft=%g dt=%g \n",nt,ft,dt);
 	fprintf(jpfp," dxm=%g fmax=%g\n",dxm,fmax);
 	fprintf(jpfp," noff=%d off0=%g doff=%g\n",noff,off0,doff);
	fprintf(jpfp," v0=%g dvz=%g \n",v0,dvz);
 	fprintf(jpfp," aperx=%g offmax=%g angmax=%g\n",aperx,offmax,angmax);
 	fprintf(jpfp," ntr=%d mtr=%d ls=%d\n",ntr,mtr,ls);
	fprintf(jpfp," ================\n");
	fflush(jpfp);

	/* allocate space */
	mig = ealloc3float(nzo,nxo,noff);
	ttab = ealloc3float(nzt,nxt,ns);
	tt = ealloc2float(nzt,nxt);
	tsum = ealloc2float(nzt,nxt);
        tbs = alloc2float(nzt,nr);
        tbg = alloc2float(nzt,nr);
        pbs = alloc2float(nzt,nr);
        pbg = alloc2float(nzt,nr);
        sigbs = alloc2float(nzt,nr);
        sigbg = alloc2float(nzt,nr);
        angbs = alloc2float(nzt,nr);
        angbg = alloc2float(nzt,nr);


 	memset((void *) mig[0][0],0,noff*nxo*nzo*sizeof(float)); 

 	fprintf(jpfp," input traveltime tables \n");
                       
	/* compute traveltime residual	*/
	for(is=0; is<ns; ++is){
		nseek = (off_t) nxt*nzt*is;
		efseeko(ttfp,nseek*((off_t) sizeof(float)),SEEK_SET);
		fread(ttab[is][0],sizeof(float),nxt*nzt,ttfp);
		s = fs+is*ds;
		resit(nxt,fxt,dxt,nzt,nr,dxo,tb[is],ttab[is],s);
		
	}

	fprintf(jpfp," start migration ... \n");
	fprintf(jpfp," \n");
	fflush(jpfp);
	

	jtr = 1;
	ktr = 0;
	do {
		/* determine offset index	*/
	    float as,res;
	    sx = tr.sx;
	    gx = tr.gx;
	    io = (int)((gx-sx-off0)/doff+0.5);
	    if(io<0) io = 0;
	    if(io>=noff) io = noff-1;

	    if(MIN(sx,gx)>=fs && MAX(sx,gx)<=es && 
	       MAX(gx-sx,sx-gx)<=offmax ){
		/*     migrate this trace	*/
	    	as = (sx-fs)/ds;
	    	is = (int)as;
		if(is==ns-1) is=ns-2;
		res = as-is;
		if(res<=0.01) res = 0.0;
		if(res>=0.99) res = 1.0;
		sum2(nxt,nzt,1-res,res,ttab[is],ttab[is+1],tsum);
                sum2(nr,nzt,1-res,res,pb[is],pb[is+1],pbs);
                sum2(nr,nzt,1-res,res,tb[is],tb[is+1],tbs);
                sum2(nr,nzt,1-res,res,sigb[is],sigb[is+1],sigbs);
                sum2(nr,nzt,1-res,res,angb[is],angb[is+1],angbs);

	    	as = (gx-fs)/ds;
	    	is = (int)as;
		if(is==ns-1) is=ns-2;
		res = as-is;
		if(res<=0.01) res = 0.0;
		if(res>=0.99) res = 1.0;
		sum2(nxt,nzt,1-res,res,ttab[is],ttab[is+1],tt);
		sum2(nxt,nzt,1,1,tt,tsum,tsum);
                sum2(nr,nzt,1-res,res,pb[is],pb[is+1],pbg);
                sum2(nr,nzt,1-res,res,tb[is],tb[is+1],tbg);
                sum2(nr,nzt,1-res,res,sigb[is],sigb[is+1],sigbg);
                sum2(nr,nzt,1-res,res,angb[is],angb[is+1],angbg);

                as = (sx-fxi)/dxi;
                is = (int)as;
                if(is==nxi-1) is=nxi-2;
                res = as-is;
                if(res<=0.01) res = 0.0;
                if(res>=0.99) res = 1.0;
                nangls = (1-res)*nangl[is] + res*nangl[is+1];

                as = (gx-fxi)/dxi;
                is = (int)as;
                if(is==nxi-1) is=nxi-2;
                res = as-is;
                if(res<=0.01) res = 0.0;
                if(res>=0.99) res = 1.0;
                nanglg = (1-res)*nangl[is] + res*nangl[is+1];

		mig2d(tr.data,nt,ft,dt,sx,gx,mig[io],aperx,
		  nxo,fxo,dxo,nzo,fzo,dzo,
		  ls,mtmax,dxm,fmax,angmax,tbs,tbg,pbs,pbg,angbs,angbg,
		  nr,tsum,nzt,fzt,dzt,nxt,fxt,dxt,szif,nangls,nanglg,
		  sigbs,sigbg);

	        ktr++;
	        if((jtr-1)%mtr ==0 ){
			fprintf(jpfp," migrated trace %d\n",jtr);
			fflush(jpfp);
	    	}
	    }
	    jtr++;
	} while (fgettr(infp,&tr) && jtr<ntr);

	fprintf(jpfp," migrated %d traces in total\n",ktr);

	memset((void *) &tro, 0, 240);
	tro.ns = nzo;
	tro.d1 = dzo;
	tro.dt = 1000*(int)(1000*dt+0.5);
	tro.delrt = 0.0;
	tro.f1 = fzo;
	tro.f2 = fxo;
	tro.d2 = dxo;
	tro.trid = 200;

	scal = 4/sqrt(PI)*dxm/v0;
	for(ixo=0; ixo<nxo; ixo++) {
		for(io=0; io<noff; io++)  {
			memcpy((void *) tro.data,
				(const void *) mig[io][ixo], nzo*sizeof(float));
			tro.offset = off0+io*doff;
			tro.tracr = 1+ixo;
			tro.tracl = 1+io+noff*ixo;
			tro.cdp = fxo+ixo*dxo;
			tro.cdpt = 1+io;

			for(izo=0; izo<nzo; ++izo)
				tro.data[izo] *=scal;

			/* write out */
			fputtr(outfp,&tro); 

		}
	}

	fprintf(jpfp," \n");
	fprintf(jpfp," output done\n");
	fflush(jpfp);

	efclose(jpfp);
	efclose(outfp);

	    
	free2float(tsum);
	free2float(tt);
        free2float(tbs);
        free2float(tbg);
        free2float(pbs);
        free2float(pbg);
        free2float(sigbs);
        free2float(sigbg);
        free2float(angbs);
        free2float(angbg);
	free3float(pb);
	free3float(tb);
	free3float(sigb);
	free3float(angb);
	free3float(ttab);
	free3float(mig);
	return(CWP_Exit());
}

/**********************************************************************
        Subroutines adapted from Dave Hale's modeling library

        decodeReflectors
        decodeReflector
        breakreflectors
        makeref                 Autor: Dave Hale, CSM, 09/17/91
**********************************************************************/

void decodeSurfaces(int *nrPtr,int **nxzPtr,float ***xPtr,float ***zPtr)
/*************************************************************************
decodeSurfaces - parse surfaces parameter string
**************************************************************************
Output:
nrPtr           pointer to nr an int specifying number of surfaces = 2
nxzPtr          pointer to nxz specifying number of (x,z) pairs defining the
                surfaces
xPtr            pointer to array[x][nr] of x values for each surface
zPtr            array[z][nr] of z values for each surface
**************************************************************************/
{
        int nr,*nxz,ir;
        float **x,**z;
	char t[6144],*s;

        /* count surfaces */
        nr = countparname("surf");
        if (nr==0) nr = 1;

        /* allocate space */
        nxz = ealloc1(nr,sizeof(int));
        x = ealloc1(nr,sizeof(float*));
        z = ealloc1(nr,sizeof(float*));

        /* get surfaces */
        for (ir=0; ir<nr; ++ir) {
                if (!getnparstring(ir+1,"surf",&s)) s = "0,0;99999,0";
                strcpy(t,s);
                if (!decodeSurface(t,&nxz[ir],&x[ir],&z[ir]))
                        err("Surface number %d specified "
                                "incorrectly!\n",ir+1);
        }

        /* set output parameters before returning */
        *nrPtr = nr;
        *nxzPtr = nxz;
        *xPtr = x;
        *zPtr = z;
}

/* parse one surface specification; return 1 if valid, 0 otherwise */

int decodeSurface (char *string,int *nxzPtr, float **xPtr, float **zPtr)
/**************************************************************************
decodeSurface - parse one surface specification
***************************************************************************
Input:
string          string representing surface

Output:
nxzPtr          pointer to number of x,z pairs
xPtr            array of x values for one surface
zPtr            array of z values for one surface
**************************************************************************/
{
        int nxz,ixz;
        float *x,*z;
        char *s,*t;

        s = string;
        s = strtok(s,",;\0");
        /* count x and z values, while splitting string into tokens */
        for (t=s,nxz=0; t!=NULL; ++nxz)
                t = strtok(NULL,",;\0");

        /* total number of values must be even */
        if (nxz%2) return 0;

        /* number of (x,z) pairs */
        nxz /= 2;

        /* 2 or more (x,z) pairs are required */
        if (nxz<2) return 0;

        /* allocate space */
        x = ealloc1(nxz,sizeof(float));
        z = ealloc1(nxz,sizeof(float));

        /* convert (x,z) values */
        for (ixz=0; ixz<nxz; ++ixz) {
                x[ixz] = atof(s);
                s += strlen(s)+1;
                z[ixz] = atof(s);
                s += strlen(s)+1;
        }

        /* set output parameters before returning */
        *nxzPtr = nxz;
        *xPtr = x;
        *zPtr = z;
        return 1;
}

void makesurf(float dsmax,int nr,int *nu,float **xu,float **zu,
        Surface **r)
/*****************************************************************************
Make piecewise cubic surfaces
******************************************************************************
Input:
dsmax           maximum length of surface segment
nr              number of surfaces = 2
nu              array[nr] of numbers of (x,z) pairs; u = 0, 1, ..., nu[ir]
xu              array[nr][nu[ir]] of surface x coordinates x(u)
zu              array[nr][nu[ir]] of surface z coordinates z(u)

Output:
r               array[nr] of surfaces
******************************************************************************
Notes:
Space for the nu, xu, and zu arrays is freed by this function, since
they are only used to construct the surfaces.
*****************************************************************************/
{
        int ir,iu,nuu,iuu,ns,is;
        float x,z,xlast,zlast,dx,dz,duu,uu,ds,fs,rsx,rsz,rsxd,rszd,
                *u,*s,(*xud)[4],(*zud)[4],*us;
        SurfaceSegment *ss;
        Surface *rr;

        /* allocate space for surfaces */
        *r = rr = ealloc1(nr,sizeof(Surface));

        /* loop over surfaces */
        for (ir=0; ir<nr; ++ir) {

                /* compute cubic spline coefficients for uniformly sampled u */
                u = ealloc1float(nu[ir]);
                for (iu=0; iu<nu[ir]; ++iu)
                        u[iu] = iu;
                xud = (float(*)[4])ealloc1float(4*nu[ir]);
                csplin(nu[ir],u,xu[ir],xud);
                zud = (float(*)[4])ealloc1float(4*nu[ir]);
                csplin(nu[ir],u,zu[ir],zud);

                /* finely sample x(u) and z(u) and compute length s(u) */
                nuu = 20*nu[ir];
                duu = (u[nu[ir]-1]-u[0])/(nuu-1);
                s = ealloc1float(nuu);
                s[0] = 0.0;
                xlast = xu[ir][0];
                zlast = zu[ir][0];
                for (iuu=1,uu=duu; iuu<nuu; ++iuu,uu+=duu) {
                        intcub(0,nu[ir],u,xud,1,&uu,&x);
                        intcub(0,nu[ir],u,zud,1,&uu,&z);
                        dx = x-xlast;
                        dz = z-zlast;
                        s[iuu] = s[iuu-1]+sqrt(dx*dx+dz*dz);
                        xlast = x;
                        zlast = z;
                }

                /* compute u(s) from s(u) */
                ns = 1+s[nuu-1]/dsmax;
                ds = s[nuu-1]/ns;
                fs = 0.5*ds;
                us = ealloc1float(ns);
                yxtoxy(nuu,duu,0.0,s,ns,ds,fs,0.0,(float)(nu[ir]-1),us);

                /* compute surface segments uniformly sampled in s */
                ss = ealloc1(ns,sizeof(SurfaceSegment));
                for (is=0; is<ns; ++is) {
                        intcub(0,nu[ir],u,xud,1,&us[is],&rsx);
                        intcub(0,nu[ir],u,zud,1,&us[is],&rsz);
                        intcub(1,nu[ir],u,xud,1,&us[is],&rsxd);
                        intcub(1,nu[ir],u,zud,1,&us[is],&rszd);
                        ss[is].x = rsx;
                        ss[is].z = rsz;
                        ss[is].c = rsxd/sqrt(rsxd*rsxd+rszd*rszd);
                        ss[is].s = -rszd/sqrt(rsxd*rsxd+rszd*rszd);
                }

                /* fill in surface structure */
                rr[ir].ns = ns;
                rr[ir].ds = ds;
                rr[ir].ss = ss;

                /* free workspace */
                free1float(us);
                free1float(s);
                free1float(u);
                free1float((float*)xud);
                free1float((float*)zud);

                /* free space replaced by surface segments */
                free1(xu[ir]);
                free1(zu[ir]);
        }

        /* free space replaced by surface segments */
        free1(nu);
        free1(xu);
        free1(zu);
}

/* Estimation of the Z coordinate of the surfaces  */
  void zcoorSurfaces(float fx,float dx,int nx,float fxs,float dxs,int nxs,
        Surface *srf,float *szif,float *sz,float *nangl)
/*****************************************************************************
From the cubic spline calculation, the Z coor. of the surface are estimated
******************************************************************************
Input:
fx              First x position in the surfaces
dx              x sampling interval
nx              number of output location points
fxs             x-coordinate of the first shot
dxs             x-coordinate increment of shots
nxs             number of shots
srf             surface structure

Output:
szif            z[] coordinates of the surface (output locations)
sz              z[] coordinates of the current surface (shot positions)
nangl           array of angles that the normal form with the vertical
*****************************************************************************/
{
        int   ik,jx,ns,ixi,ix;
        float x,ax,ax0,az,temp;
        float *ssx,*ssz,*sss;
        SurfaceSegment *ss;

        ns = srf[0].ns;
        ssx = alloc1float(ns);
        ssz = alloc1float(ns);
        sss = alloc1float(ns);

        /* loop over surface  */
        /* number of segments, segment length */
        ss = srf[0].ss;
        for (jx=0; jx<ns ; jx++) {
                ssx[jx] = ss[jx].x;
                ssz[jx] = ss[jx].z;
                sss[jx] = ss[jx].s;
        }
        ixi = 0;
        for (ix=0; ix<nx; ++ix) {
            x = fx + ix*dx;
            for (ik=ixi; ik<ns-1; ++ik) {
                if (ssx[ik]<=x && ssx[ik+1]>=x) {
                   az = ssz[ik] - ssz[ik+1];
                   ax0 = ssx[ik+1] - x;
                   ax = ssx[ik+1] - ssx[ik];
                   szif[ix] = ax0*az/ax+ssz[ik+1];
                   az = sss[ik] - sss[ik+1];
                   temp = ax0*az/ax + sss[ik+1];
                   nangl[ix] = asin(temp);
                   ixi = ik;
		}
            }
        }
        ixi = 0;
        for (ix=0; ix<nxs; ++ix) {
            x = fxs + ix*dxs;
            for (ik=ixi; ik<ns-1; ++ik) {
                if (ssx[ik]<=x && ssx[ik+1]>=x) {
                   az = ssz[ik] - ssz[ik+1];
                   ax0 = ssx[ik+1] - x;
                   ax = ssx[ik+1] - ssx[ik];
                   sz[ix] = ax0*az/ax+ssz[ik+1];
		   ixi = ik;
                }
            }
        }
        szif[0] = 2*szif[1] - szif[2];
        szif[nx-1] = 2*szif[nx-2] - szif[nx-3];
        sz[0] = 2*sz[1] - sz[2];
        nangl[0] = 2*nangl[1] - nangl[2];
        sz[nxs-1] = 2*sz[nxs-2] - sz[nxs-3];
        nangl[nx-1] = 2*nangl[nx-2] - nangl[nx-3];

        free1(ssx);
        free1(ssz);
        free1(sss);
}

/* residual traveltime calculation based  on reference   time	*/
  void resit(int nx,float fx,float dx,int nz,int nr,float dr,
		float **tb,float **t,float x0)
{
	int ix,iz,jr;
	float xi,ar,sr,sr0;

	for(ix=0; ix<nx; ++ix){
		xi = fx+ix*dx-x0;
		ar = abs(xi)/dr;
		jr = (int)ar;
		sr = ar-jr;
		sr0 = 1.0-sr;
		if(jr>nr-2) jr = nr-2;
		for(iz=0; iz<nz; ++iz)
			t[ix][iz] -= sr0*tb[jr][iz]+sr*tb[jr+1][iz];
	}
} 

/* lateral interpolation	*/

/* sum of two tables	*/
  void sum2(int nx,int nz,float a1,float a2,float **t1,float **t2,float **t)
{
	int ix,iz;

	for(ix=0; ix<nx; ++ix) 
		for(iz=0; iz<nz; ++iz)
			t[ix][iz] = a1*t1[ix][iz]+a2*t2[ix][iz];
}

/* compute  reference traveltime and slowness   */
      void timeb(int nr,int nz,float dr,float dz,float fz,float z0,float a,
        float v0,float **t,float **p,float **sig,float **ang)
{
        int  ir,iz;
        float r,z,v,rc,oa,temp,rou,zc;


        if( a==0.0) {
                for(ir=0,r=0;ir<nr;++ir,r+=dr)
                        for(iz=0,z=fz-z0;iz<nz;++iz,z+=dz){
                                rou = sqrt(r*r+z*z);
                                if(rou<dz) rou = dz;
                                t[ir][iz] = rou/v0;
                                p[ir][iz] = r/(rou*v0);
                                sig[ir][iz] = v0*rou;
                                ang[ir][iz] = asin(r/rou);
                        }
        } else {
                oa = 1.0/a;     zc = v0*oa;
                for(ir=0,r=0;ir<nr;++ir,r+=dr)
                        for(iz=0,z=fz+zc-z0;iz<nz;++iz,z+=dz){
                                rou = sqrt(r*r+z*z);
                                v = v0+a*(z-zc+z0);
                                if(ir==0){
                                        t[ir][iz] = log(v/v0)*oa;
                                        p[ir][iz] = 0.0;
                                        ang[ir][iz] = 0.0;
                                        sig[ir][iz] = 0.5*(z+0.1*dz-zc+z0)
						*(v0+v);
                                } else {
                                        rc = (r*r+z*z-zc*zc)/(2.0*r);
                                        rou = sqrt(zc*zc+rc*rc);
                                        t[ir][iz] = log((v*(rou+rc))
                                                /(v0*(rou+rc-r)))*oa;
                                        p[ir][iz] = sqrt(rou*rou-rc*rc)
                                                /(rou*v0);
                                        temp = v0*p[ir][iz];
                                        if(temp>1.0) temp = 1.0;
                                        ang[ir][iz] = asin(temp);
                                        sig[ir][iz] = a*rou*r;
                                }
                        }
        }
}

void filt(float *trace,int nt,float dt,float fmax,int ls,int m,float *trf);

  void mig2d(float *trace,int nt,float ft,float dt,
	float sx,float gx,float **mig,float aperx,
  	int nx,float fx,float dx,float nz,float fz,float dz,
	int ls,int mtmax,float dxm,float fmax,float angmax,
        float **tbs,float **tbg,float **pbs,
        float **pbg,float **angbs,float **angbg,int nr,float **tsum,
        int nzt,float fzt,float dzt,int nxt,float fxt,float dxt,float *szif,
        float nangls,float nanglg,float **sigbs,float **sigbg)
/*****************************************************************************
Migrate one trace 
******************************************************************************
Input:
*trace		one seismic trace 
nt		number of time samples in seismic trace
ft		first time sample of seismic trace
dt		time sampleing interval in seismic trace
sx,gx		lateral coordinates of source and geophone 
aperx		lateral aperature in migration
nx,fx,dx,nz,fz,dz  dimension parameters of migration region
ls		=1 for line source; =0 for point source
mtmax		number of time samples in triangle filter
dxm		midpoint sampling interval
fmax		frequency-highcut for input trace	 
angmax		migration angle aperature from vertical 	 
(tb,pb,sigb,angb)s,g	reference traveltime, lateral slowness, cosine of 
		incident angle, and emergent angle
nr		number of lateral samples in reference quantities
tsum		sum of residual traveltimes from shot and receiver
nxt,fxt,dxt,nzt,fzt,dzt		dimension parameters of traveltime table
szif            array[] of the z component of the recording surface (ignored)
nangls,nanglg   Angles that the normal unit vector of the current surface
                forms with the vertical at source and geophone locations.
Output:
mig		migrated section
*****************************************************************************/
{
	int nxf,nxe,nxtf,nxte,ix,iz,iz0,izt0,nzp,jrs,jrg,jz,jt,mt,jx;
	float xm,x,dis,rxz,ar,srs,srg,srs0,srg0,sigp,z0,rdz,ampd,res0,
	      angs,angg,ax,ax0,pmin,odt=1.0/dt,pd,az,sz,sz0,sigmas,sigmag,
	      at,td,res,temp,sigs,sigg;
	float **tmt,**ampt,**ampti,*tm,*amp,*ampi,*tzt,*trf,*zpt;

	tmt = alloc2float(nzt,nxt);
	ampt = alloc2float(nzt,nxt);
	ampti = alloc2float(nzt,nxt);
	tm = szif; /* dummy to keep compiler happy */
	tm = alloc1float(nzt);
	tzt = alloc1float(nzt);
	amp = alloc1float(nzt);
	ampi = alloc1float(nzt);
	zpt = alloc1float(nxt);
	trf = alloc1float(nt+2*mtmax);

	z0 = (fz-fzt)/dzt;
	rdz = dz/dzt;
	pmin = 1.0/(2.0*dxm*fmax);
	
	filt(trace,nt,dt,fmax,ls,mtmax,trf);

	xm = 0.5*(sx+gx);
	rxz = (angmax==90)?0.0:1.0/tan(angmax*PI/180.);
	nxtf = (xm-aperx-fxt)/dxt;
	if(nxtf<0) nxtf = 0;
	nxte = (xm+aperx-fxt)/dxt+1;
	if(nxte>=nxt) nxte = nxt-1;

	/* compute amplitudes and filter length	*/
	for(ix=nxtf; ix<=nxte; ++ix){
		x = fxt+ix*dxt;
		dis = (xm>=x)?xm-x:x-xm;
		izt0 = ((dis-dxt)*rxz-fzt)/dzt-1;
		if(izt0<0) izt0 = 0;
		if(izt0>=nzt) izt0 = nzt-1;

		ar = (sx>=x)?(sx-x)/dx:(x-sx)/dx;
		jrs = (int)ar;
		if(jrs>nr-2) jrs = nr-2;
		srs = ar-jrs;
		srs0 = 1.0-srs;
		ar = (gx>=x)?(gx-x)/dx:(x-gx)/dx;
		jrg = (int)ar;
		if(jrg>nr-2) jrg = nr-2;
		srg = ar-jrg;
		srg0 = 1.0-srg;
		sigp = ((sx-x)*(gx-x)>0)?1.0:-1.0;
                sigs = ((sx-x)<0)?1.0:-1.0;
                sigg = ((gx-x)<0)?1.0:-1.0;
		zpt[ix] = fzt+(nzt-1)*dzt;

		for(iz=izt0; iz<nzt; ++iz){
			angs = srs0*angbs[jrs][iz]+srs*angbs[jrs+1][iz];
                        angs = ABS(sigs*angs - nangls);
			angg = srg0*angbg[jrg][iz]+srg*angbg[jrg+1][iz];
                        angg = ABS(sigg*angg - nanglg);
			sigmas = srs0*sigbs[jrs][iz]+srs*sigbs[jrs+1][iz]; 
			sigmag = srg0*sigbg[jrg][iz]+srg*sigbg[jrg+1][iz]; 
                        ampd = cos(angs)*sigmag + cos(angg)*sigmas;
			ampd /= (cos(nangls)*cos(nanglg)*sqrt(sigmas*sigmag));

        /* Filter of 90 degrees to the operator (Peels,1988)

           Peels, G. L., 1988, True amplitude wave field extrapolation
           with applications in  seismic shot record redatuming, PhD
           thesis, Delft University of Technology.              */

                        if(ABS(angs)>=(PI/2) || ABS(angg)>=(PI/2)) {
                                ampd=0.0;
                        }

			if(ampd<0.0) ampd = -ampd;
			ampt[ix][iz] = ampd;

			pd = srs0*pbs[jrs][iz]+srs*pbs[jrs+1][iz]+sigp 
			     *(srg0*pbg[jrg][iz]+srg*pbg[jrg+1][iz]);
			if(pd<0.0) pd = -pd;
			temp = pd*dxm*odt;
			if(temp<1) temp = 1.0;
			if(temp>mtmax) temp = mtmax;
			ampti[ix][iz] = ampd/(temp*temp);
			tmt[ix][iz] = temp;
			if(pd<pmin && zpt[ix]>fzt+(nzt-1.1)*dzt) 
				zpt[ix] = fzt+iz*dzt;

		}
	}

	nxf = (xm-aperx-fx)/dx+0.5;
	if(nxf<0) nxf = 0;
	nxe = (xm+aperx-fx)/dx+0.5;
	if(nxe>=nx) nxe = nx-1;
	
	/* interpolate amplitudes and filter length along lateral	*/
	for(ix=nxf; ix<=nxe; ++ix){
		x = fx+ix*dx;
		dis = (xm>=x)?xm-x:x-xm;
		izt0 = (dis*rxz-fzt)/dzt;
		if(izt0<0) izt0 = 0;
		if(izt0>=nzt) izt0 = nzt-1;
		iz0 = (dis*rxz-fz)/dz;
		if(iz0<0) iz0 = 0;
		if(iz0>=nz) iz0 = nz-1;

		ax = (x-fxt)/dxt;
		jx = (int)ax;
		ax = ax-jx;
		if(ax<=0.01) ax = 0.;
		if(ax>=0.99) ax = 1.0;
		ax0 = 1.0-ax;
		if(jx>nxte-1) jx = nxte-1;
		if(jx<nxtf) jx = nxtf;

		ar = (sx>=x)?(sx-x)/dx:(x-sx)/dx;
		jrs = (int)ar;
		if(jrs>nr-2) jrs = nr-2;
		srs = ar-jrs;
		srs0 = 1.0-srs;
		ar = (gx>=x)?(gx-x)/dx:(x-gx)/dx;
		jrg = (int)ar;
		if(jrg>nr-2) jrg = nr-2;
		srg = ar-jrg;
		srg0 = 1.0-srg;

		for(iz=izt0; iz<nzt; ++iz){
		    tzt[iz] = ax0*tsum[jx][iz]+ax*tsum[jx+1][iz]
				+srs0*tbs[jrs][iz]+srs*tbs[jrs+1][iz]
				+srg0*tbg[jrg][iz]+srg*tbg[jrg+1][iz];

		    amp[iz] = ax0*ampt[jx][iz]+ax*ampt[jx+1][iz];
		    ampi[iz] = ax0*ampti[jx][iz]+ax*ampti[jx+1][iz];
		    tm[iz] = ax0*tmt[jx][iz]+ax*tmt[jx+1][iz];

		}

		nzp = (ax0*zpt[jx]+ax*zpt[jx+1]-fz)/dz+1.5;
		if(nzp<iz0) nzp = iz0;
		if(nzp>nz) nzp = nz;

		/* interpolate along depth if operater aliasing 	*/
		for(iz=iz0; iz<nzp; ++iz) {
			az = z0+iz*rdz;
			jz = (int)az;
			if(jz>=nzt-1) jz = nzt-2;
			sz = az-jz;
			sz0 = 1.0-sz;
			td = sz0*tzt[jz]+sz*tzt[jz+1];
			at = (td-ft)*odt+mtmax;
			jt = (int)at;
			if(jt > mtmax && jt < nt+mtmax-1){
			    ampd = sz0*ampi[jz]+sz*ampi[jz+1];
			    mt = (int)(0.5+sz0*tm[jz]+sz*tm[jz+1]);
			    res = at-jt;
			    res0 = 1.0-res;
 			    temp = (res0*(-trf[jt-mt]+2.0*trf[jt]-trf[jt+mt]) 
				+res*(-trf[jt-mt+1]+2.0*trf[jt+1]
				-trf[jt+mt+1]))*ampd;				
			    mig[ix][iz] += temp;

			}
		}

		/* interpolate along depth if not operater aliasing 	*/
		for(iz=nzp; iz<nz; ++iz) {
			az = z0+iz*rdz;
			jz = (int)az;
			if(jz>=nzt-1) jz = nzt-2;
			sz = az-jz;
			sz0 = 1.0-sz;
			td = sz0*tzt[jz]+sz*tzt[jz+1];
			at = (td-ft)*odt;
			jt = (int)at;
			if(jt > 0 && jt < nt-1){
			    ampd = sz0*amp[jz]+sz*amp[jz+1];
			    res = at-jt;
			    res0 = 1.0-res;
 			    temp = (res0*trace[jt]+res*trace[jt+1])*ampd; 
			    mig[ix][iz] += temp;
			}
		}

	}

	free2float(ampt);
	free2float(ampti);
	free2float(tmt);
	free1float(amp);
	free1float(ampi);
	free1float(zpt);
	free1float(tm);
	free1float(tzt);
	free1float(trf);
}

void filt(float *trace,int nt,float dt,float fmax,int ls,int m,float *trf)
/*********************************************************************
  Low-pass filter, integration and phase shift for input data	 
  input: 
    trace(nt)	single seismic trace
    fmax	high cut frequency
    ls		ls=1, line source; ls=0, point source
  output:
    trace(nt) 	filtered and phase-shifted seismic trace 
    tracei(nt) 	filtered, integrated and phase-shifted seismic trace 
**********************************************************************/
{
	static int nfft=0, itaper, nw, nwf;
	static float *taper, *amp, *ampi, dw;
	int  it, iw, itemp;
	float temp, ftaper, const2, *rt;
	complex *ct;

	fmax *= 2.0*PI;
	ftaper = 0.1*fmax;
	const2 = 0.5*sqrt(2.0);

	if(nfft==0) {
        	/* Set up FFT parameters */
        	nfft = npfaro(nt+2*m, 2*(nt+2*m));
        	if (nfft >= SU_NFLTS || nfft >= 720720)
                	err("Padded nt=%d -- too big", nfft);

        	nw = nfft/2 + 1;
		dw = 2.0*PI/(nfft*dt);

		itaper = 0.5+ftaper/dw;
		taper = ealloc1float(2*itaper+1);
		for(iw=-itaper; iw<=itaper; ++iw){
			temp = (float)iw/(1.0+itaper); 
			taper[iw+itaper] = (1-temp)*(1-temp)*(temp+2)/4;
		}

		nwf = 0.5+fmax/dw;
		if(nwf>nw-itaper-1) nwf = nw-itaper-1;
		amp = ealloc1float(nwf+itaper+1);
		ampi = ealloc1float(nwf+itaper+1);
		amp[0] = ampi[0] = 0.;
		for(iw=1; iw<=nwf+itaper; ++iw){
			amp[iw] = sqrt(dw*iw)/nfft;
			ampi[iw] = 0.5/(1-cos(iw*dw*dt));
		}
	}

        /* Allocate fft arrays */
        rt   = ealloc1float(nfft);
        ct   = ealloc1complex(nw);

        memcpy(rt, trace, nt*FSIZE);
        memset((void *) (rt + nt), 0, (nfft-nt)*FSIZE); 
        pfarc(1, nfft, rt, ct);

	for(iw=nwf-itaper;iw<=nwf+itaper;++iw){
		itemp = iw-(nwf-itaper);
		ct[iw].r = taper[itemp]*ct[iw].r; 
		ct[iw].i = taper[itemp]*ct[iw].i; 
	}
	for(iw=nwf+itaper+1;iw<nw;++iw){
		ct[iw].r = 0.; 
		ct[iw].i = 0.; 
	}

       	if(!ls){
		for(iw=0; iw<=nwf+itaper; ++iw){
			/* phase shifts PI/4 + PI/4 (Half derivative)  */
			temp = -ct[iw].i*amp[iw];
			ct[iw].i = ct[iw].r*amp[iw];
			ct[iw].r = temp;
		    }
	} else {
		for(iw=0; iw<=nwf+itaper; ++iw){
                        /* phase shifts PI/4 -Half derivative   */
                        temp = (ct[iw].r-ct[iw].i)*amp[iw]*const2;
                        ct[iw].i = (ct[iw].r+ct[iw].i)*amp[iw]*const2;
                        ct[iw].r = temp;
		}
	}              
        pfacr(-1, nfft, ct, rt);
		
        /* Load traces back in */
	for (it=0; it<nt; ++it) trace[it] = rt[it];

        /* Integrate traces   */
	for(iw=0; iw<=nwf+itaper; ++iw){
		ct[iw].i = ct[iw].i*ampi[iw];
		ct[iw].r = ct[iw].r*ampi[iw];
	}
        pfacr(-1, nfft, ct, rt);
        for (it=0; it<m; ++it)  trf[it] = rt[nfft-m+it];
        for (it=0; it<nt+m; ++it)  trf[it+m] = rt[it];

	free1float(rt);
	free1complex(ct);
}
